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We demonstrate that Monte Carlo sampling can be used to efficiently extract the expectation value 
of projected entangled pair states with a large virtual bond dimension. We use the simple update rule 
introduced by Xiang et al. in Phys. Rev. Lett 101, 090603 (2008) to obtain the tensors describing 
the ground state wavefunction of the Antiferromagnetic Heisenberg model and evaluate the finite 
size energy and staggered magnetization for square lattices with periodic boundary conditions of 
linear sizes up to L — 16 and virtual bond dimensions up to D = 16. The finite size magnetization 
errors are 0.003(2) and 0.013(2) at D = 16 for a system of size L = 8, 16 respectively. Finite D 
extrapolation provides exact finite size magnetization for L = 8, and reduces the magnetization 
error to 0.005(3) for L = 16, significantly improving the previous state of the art results. 

PACS numbers: 02.70.Ss, 75.10.Jm, 75.40.Mg, 75.40.Cx 
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The efficient simulation of strongly correlated quan- 
tum many body systems presents one of the major open 
problems and challenges in condensed matter physics. A 
major step forward was made by Steven White — in the 
case of 1 dimensional quantum spin chains by introduc- 
ing the density matrix renormalization group (DMRG), 
which soon became the method of choice for simulating 
1 dimensional manybody systems at zero temperature. 
By reformulating DMRG as a variational method within 
the class of matrix product states (MPS)2r— , it has be- 
come clear how DMRG can be generalized to deal with 
systems in two dimensions^; the quantum states of the 
corresponding variational class are known as projected 
entangled pair states (PEPS) and are part of the class 
called tensor product states which also includes the mul- 
tiscale entanglement renormalization ansatz 8 and infinite 
PEPS 9 . More recently, it has also been demonstrated 
how the PEPS class can take into account fermionic 
anti-commutation relations — . Numerical algorithms 
based on these ansatze, such as variational minimization 
of the ground state energy and imaginary time evolution 
are also developing fast — il & 17 i 18 , and a wide range of 
applications has been studied 1 ^— . 

The computational complexity of algorithms based on 
the PEPS ansatz with virtual bond dimension D scales as 
D 12 for the finite PEPS algorithm with open boundary 
condition 1 !, x 3 D 4 for the infinite PEPS (iPEPS) algo- 
rithm*, x e f° r the tensor entanglement renormalization 
(TERG) algorithm for square lattices^ and x 5 for hon- 
eycomb lattice s 1 ' 18 , where x is the number of Schmidt co- 
efficients kept in the various approximations. The large 
scaling power presents the main bottleneck in scaling up 
the number of variational parameters, which is neces- 
sary near second order phase transitions — . The com- 
mon characteristics of all these algorithms is that the 
tensor network is always contracted over the physical in- 
dices, which effectively squares the computational cost 
of contracting the tensor network as compared to a ten- 
sor network corresponding to a classical spin system. As 
first shown i n 29 ' 30 for the case of matrix product states 
and string bond states, a square root speed up can be 



obtained by using importance sampling over the physi- 
cal indices. We will show how to adapt an importance 
sampling technique to PEPS. The efficiency will depend 
on the contraction algorithm chosen. In this paper we 
demonstrate it using the TERG method. 

The Antiferromagnetic Heisenberg model on a square 
lattice with length L has been well studied by stochas- 
tic series expansions (SSE)21, however it is notori- 
ously hard for tensor network wavefunctions to precisely 
capture the ground state order parameter (staggered 
magnetization) 25 . Various attempts have been made to 
extract the right magnetization, e.g. using iPEPS algo- 
rithm on square lattice 2 ^ and the second renormalization 
of tensor network state (SRG) on honeycomb lattice 2 ^. 
However, all those attempts indicated that a tensor prod- 
uct state (TPS) with a finite D has much larger stag- 
gered magnetization in the thermodynamic limit. The 
main reason for that is probably the fact that all TPS 
methods favour states with a small amount of entangle- 
ment, and a larger local order parameter indeed leads to 
states with a smaller amount of entanglement due to the 
monogamy property of entanglement 34 . 

The simple update proposed in RefJ« is an extremely 
fast imaginary time evolution (projection) method, which 
makes a simple estimation of the entanglement between 
the sub-system and the environment and integrates it 
in the evolution step. The evolution does not aim at 
the time dependent state at imaginary time t; it aims 
at that, in the long run, the accumulative effect of many 
non-sufficient improvements will eventually drive the sys- 
tem to the ground state. Since there is no notion of the 
lattice size in this update, one can claim the ground state 
obtained must be that of an infinite lattice. Given a ten- 
sor product description of the wavefunction with virtual 
bond dimension D and its correlation length £(D), no 
local observable will have any notion the lattice size if 
L > £,(D). We take the tensors obtained from the sim- 
ple update describing the ground state of Antiferromag- 
netic on an infinite lattice and evaluate the finite size 
energy and staggered magnetization with Monte Carlo 
(MC) sampling technique. We show that the magnetiza- 
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FIG. 1: (a) The tensor network wavefunction of a spin system 
on a square lattice, (b) The contraction of physical indices for 
calculating expectation values of a tensor network wavefunc- 
tion, (c) this results in a tensor network of bond dimension 
D 2 . 



tion indeed reaches the correct value when larger bond 
dimensions are used. 

The paper is organized as following: in Sec. U we give a 
brief introduction to the TERG algorithm, in Sec. [TTI we 
illustrate the sampling procedure using the TERG con- 
traction method, in Sec. IIIII we apply the ground state 
tensor obtained via the simple update (poorman's up- 
date) 1 - to finite size lattices and evaluate finite size expec- 
tation values via MC sampling, and finally a summary is 
given in Sec. [TV] 



I. TENSOR ENTANGLEMENT 
RENORMALIZATION ALGORITHM 

The tensor network ansatz describes quantum many- 
body states in an exponentially large Hilbert space in 
terms of local tensors T describing local degrees of free- 
dom. A graphical representation of a tensor network 
state for a spin model on a square lattice is presented 
in Fig. [Ha), for which the wavefunction is written as 



E 



tTr { T [i],si T [2], S 2 . . . T M,»w}| aii aa , . . . , s N ), 



(1) 

where TW ,Si denotes the tensor of spin s, on site i and 
|<r) = \s\, S2, ■ ■ ■ , sn) represents manybody spin configu- 
rations. The notation tTr is used to represent tensorial 
trace, generalization of the matrix trace to tensor net- 
works where tensors are traced over the virtual modes. 

If one were to calculate an expectation value for a 
given observable given a PEPS state, it would require 
contraction of a double layer tensor network obtained by 
contracting over the physical bonds first as depicted on 
Fig. |ljb) . This results in a tensor network with a squared 
bond dimension as shown on Fig. HJc). The effort for 
contracting a tensor network either as Fig. [IJa) for single 
layer or as Fig. QJc) for double layer grows exponentially 
with the system size. One way of exact contraction is to 
successively renormalize a 2 x 2 block of sites into one 
super-site, as pointed out in Ref^. Without approxi- 
mation, the dimension of the virtual bond of a super-site 




FIG. 2: (a) First decompose each 4-index tensor (in open dia- 
monds on dash lines) into two 3- index tensors (in black dots), 
then contract every 4 tensors in the shaded area into a 4-index 
tensor (in open square in (b)). (b) Repeat the decomposition- 
contraction procedure on a reduced and rotated lattice (in 
gray solid line) . (c) The total effect is that each 2x2 cluster 
on a fine lattice (in dash lines) is coarse grained into a su- 
per site (in open diamond) on a coarse grained lattice; note 
that the lattice orientation can be restored after every two 
iterations. 



in a double layer picture will grow as D 4 ™'', where n r is 
the depth of the renormalization iteration (the linear size 
of the system is L ~ 2" r ) , thus approximate contraction 
becomes necessary. One of the ways to contract the ten- 
sor network approximately is the tensor renormalization 
method 33 , which was first proposed to contract a classi- 
cal tensor network. Later on this method was general- 
ized to deal with quantum systemsi^. The contraction 
method on a square lattice can be described in Fig. [2] 
First each 4-index T-tensor is decomposed into two 3- 
index S-tensors, 
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(2) 
(3) 



where T B,A denote two ways of tensor decomposition 
according to Eq. Q © respectively. In the next step the 
four S-tensors in the shaded area in Fig. [2a) are con- 
tracted to form a coarse tensor on a reduced and rotated 
lattice, 



Eq2 c3 c4 q1 
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(4) 



This decomposition-contraction procedure can be applied 
once again on the rotated lattice (Fig. HJb)) to obtain a 
coarse lattice of half the length (Fig[2jc)), and whose 
orientation of the lattice is equal to the original one. 

A singular value decomposition (SVD) is then done to 
decompose a T-tensor into two S-tensors, 



D 



(5) 



a=l 



To prevent an exponential increase of the computational 
cost, one only keeps the largest D cut (also referred to 
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as x) singular values; this approximation maximizes the 
2-norm of vectorized T for a fixed D cut , 



Tjjkl ~ J]] Uij a K a Vklai 



(6) 



a=l 



where M denotes taking the leading D cut columns of a 
matrix M (or leading singular values of a diagonal ma- 
trix). A common strategy is to absorb the diagonal ma- 
trix A into isometries U and V to obtain the S-tensors, 
S 1 = UVX and S 3 = VVX, the same applies to S 2 , S 4 . 



II. VARIATIONAL QUANTUM MONTE CARLO 
SAMPLING AND UPDATE 



A variational quantum Monte Carlo (vQMC) method 
with tensor network states can now be based on the 
TERG contraction method to calculate the importance 
weight and the energy derivative. Following the no- 
tation of Rcf. 2 -, we extract several key equations re- 
garding measuring and updating. For a chosen con- 
figuration \a) = |si, S2» * " * t s n) we define a coefficient 
W(a) — ((j\4>), which is calculated by contracting a sin- 
gle layer tensor network as 

W(a) = tTr{T [11 ' Sl T [2] ' S2 • • • T [A,I - SJV }. (7) 

The energy expectation value reads 



where 
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E(a) 



v W{a' 



~f W(a) 
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(8) 



(9) 



The energy derivatives with respect to tensor elements 
■* are 

dE 



Tijki are obtained via 
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2<Aj$ (a)E(a))-2(A^ (a)) (E(a)), (10) 
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the () denotes Monte Carlo average and 

A^; -. . (ii) 



Ui) _ 1 dW{a) 



Let us define B(m) as the contraction of the tensor net- 
work for all sites except a site m: 



B(m) = tTrf • • T^-\Urn-x T [m+\] 



in terms of which we can express the derivative of the 
weight ((7]) with respect to Pj^i as 

dW{a) 



<9T, 



t(4) 



= ^2B(m) ijkl 6 SmMi) , 



(13) 



where we assume translation invariance symmetry, i.e. 
TW- s = T^- s (s =t,|) for all sites i,j in the lattice. 

The program starts by randomly generating a spin 
configuration \a) satisfying J^. Si = 0, i.e. we initial- 
ize our state to live in total spin sector. Given |cr), 
one initializes and stores all the intermediate T q ' p -tensor 
at the site q of the p th coarse grained lattice. Dur- 
ing the contraction, we also calculate and store scalars 
/i>p = max{|T^P|} for all T q - P , then divide T?.' p by 
/ q ' p to avoid too large or too small singular values in the 
next iteration. If we define the tensor trace of the final 
contraction step as g = tTr{T[ 1 ]<™-T[ 2 ]'™-T[ 3 ]'™''T[ 4 ]' 11 '-}, 
where n r is the number of iterations of a tensor network 
contraction (L = 2^r +1 ), then weight ([7]) can be written 
as 



W{a)=g\{r 



(14) 



q.p 



Since we do not need to update the variational param- 
eters tJ.^i in this work, we will discuss how to update 
a tensor network using energy derivatives with MC sam- 
pling technique elsewhere 4 ^. 

While describing the sampling procedure, we take the 
nearest neighbor Antiferromagnetic Heisenberg interac- 
tion as an example. Generalization to other Hamiltonian 
is straight forward. Starting from site 1 of the original 
tensor network, one looks for a pair of nearest neigh- 
bor spins that align anti-parallel with each other and flip 
them. 

The trial configuration, which we denote as |ct'), is 
accepted with probability 



P 



1, 



W 2 (a') 



W 2 {a) 



where the ratio is given by 

W(a') g' 



W(a) 



n 



f'q.p 



9 " / q 
q-p 



(15) 



(16) 
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To calculate the ratio (fTBj). one needs to recompute some 
T' q ' p tensors together with the corresponding scalars 
/' q,p and g' , store them in separate arrays for later up- 
dates. If a random number r drawn from an uniform 
distribution on the interval [0,1) satisfies r < P, the 
trial state \a'} is accepted, in which case T q,p , / q,p and 
g are replaced by T' q,p , /' q ' p and g', otherwise \a > ) is 
rejected and the original configuration \a) is kept. Mov- 
ing through all the sites on the original lattice, one at- 
tempts to flip all encountered anti-parallel pairs, accept- 
ing or rejecting according to probability (|15p. This pro- 
cedure is called a MC sweep. After each MC sweep, the 
energy and other observables of interest are measured. 
Flipping two neighboring spins does not require recom- 
puting many T' q,p tensors, which makes the contraction 
fast. On the other hand, the update is local. To reduce 
the auto correlation length, one needs to complete a MC 
sweep before making a measurement, and the computa- 
tional effort scales linearly with the system size N = L 2 . 
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FIG. 3: (a) A evolution operator is applied to the nearest 
neighbor sites A,B with weight \/Ai putting on the open 
bonds, (b) QR and LQ decompositions are done separately 
on T A and T B prior to the SVD step taken over the evolved 
bond between sites A,B. The computational cost is thus re- 
duced from D 9 to D 5 . 



III. THE ANTIFERROMAGNETIC 
HEISENBERG MODEL ON SQUARE LATTICE 

We use the simple update method of Xiang et alii to 
obtain the converged wavefunction with various virtual 
bond dimension (D = 3, 4, • • • , 20). The simple update is 
an imaginary time evolution method to obtain the ground 
state wavefunction of an infinite lattice. 

To implement the imaginary time evolution, we first 
make the Trotter decomposition of the partition function 
e~P H = [e~ T t^i H i] M for r = /3/M, then we apply an evo- 
lution operator e~ rHj to the two nearest neighbor sites 
A, B as in Fig. G2[a). It is crucial to put the weight \f~Ki 
to the open auxiliary modes to take into account the en- 
tanglement of the sub-system with the environment. The 
weights y/Xl are the singular values obtained in the previ- 
ous evolution step on the corresponding bond. According 
to Refi an SVD is done to the joint tensor T A T B to in- 
troduce a cut over the enlarged bond, and only the lead- 
ing D singular values and the corresponding left and right 
eigenvectors are kept as a projection to the sub-manifold 
of the Hilbert space where the wavefunction manifests. 
Here we made a crucial modification that drastically re- 
duces the computational cost of decomposing the joint 
matrix T A T B that is of size dD 3 x dD 3 . As illustrated in 
Fig. [3jb) , we first make QR(LQ) factorization of tensor 
T A < B ) as T A = Q A R A and T B = L B Q B , where R A and 
L B are right- and left-triangular matrices. Instead of a 
large tensor T A T B , the singular value decomposition is 
done on R A L B , which is essentially of size dD x dD, as 
R A L B = UAV T . The leading computational cost is now 
the QR(LQ) decomposition that scales only to D 5 . To 
obtain the evolved tensors on sites A and B, one has to 
remove the weight y/Xl from the decomposed tensors as 
described in Refi 

In the thermodynamic limit, the ground state of 
the Antifcrromagnetic Heisenberg model spontaneously 
breaks the SU(2) symmetry, i.e. the magnetization is 
locked in one direction. To achieve a significant accep- 
tance ratio for the Markov process with local spin flips, 
we intentionally break the SU(2) symmetry into the XY 
plane. To do this, we first attach a small (h a = 0.001J) 
staggered magnetic field in the x direction to the isotropic 



Heisenberg Hamiltonian, 
H = j£ Si ■ Sj + ha^-l^Sf, (J > 0), (17) 

here i x ,i y are x,y coordinates of site i. We update the 
tensor network wavefunction using this modified Hamil- 
tonian until it converges. Then we use the converged 
wavefunction to initialize a new update without the field. 
The Trotter steps of the imaginary time evolution is grad- 
ually reduced from t = 10 ~ 2 to 10~ 5 . A convergence 

is reached when ^^^^^ < 10" 7 , where 

T A,B (r) is the vectorized tensor at time slice r, and is 
rescaled such that the largest magnitude of the tensor 
elements is 1. We then take these converged tensors of 
various bond dimension D for the infinite lattice to com- 
pute expectation values of finite lattices with periodic 
boundary condition using MC sampling method. 

One way to define the staggered magnetization is 
through the spin-spin correlation at the longest dis- 
tance^ 

M 2 = Y,C a (L/2,L/2), (18) 

a 

where a — x,y, z, and 

C a (L/U/2) = i^S a (i I , !s )S»( !l + ^ J + ^. 

i 

(19) 

In Fig. SI we present the staggered magnetization as a 
function of inverse virtual bond dimension D for system 
sizes L = 4, 8, 16. The solid lines represents the mag- 
netization results for finite lattices obtained via SSE 31 
and resonating valence bond (RVB) projection 3 ^ meth- 
ods. For a small size L = 4, large bond dimension D > 8 
gives the exact magnetization within the statistical er- 
ror. For larger sizes L = 8, 16, the magnetization error 
at D = 16 is 0.003(2) 0.013(2) respectively. Finite D ex- 
trapolation gives the exact finite size magnetization for 
L — 8, and reduces the magnetization error to 0.005(3) 
for L = 16. 

In Fig. [5] we present the absolute error of the finite size 
energy divided by the number of bonds as a function of 
virtual bond dimension D for system sizes L = 4,8, 16. 
For all system sizes the energy error drops significantly 
at D — 10 and at D e [10 : 16] plateaus seem to set in. 

In Fig. |6] we show all three components of the spin- 
spin correlation at the longest distance C a (L/2, L/2), 
a = x,y,z. One can see that the x, y components for 
different system sizes almost fall on top of each other. 
The x, y components slightly drop at D = 5, 10 then 
followed by plateaus. The z component, on the other 
hand, largely deviates from the x and y components. For 
L = 4, the SU(2) symmetry is gradually restored with 
increasing bond dimension D; for L = 8, there is a par- 
tial growth of C z (L/2, L/2) for increasing D; however for 
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FIG. 4: Staggered magnetization as a function of 1/D for 
L — 4, 8, 16. The solid lines are finite size expectation value 
from SSE and RVB projection methods. The dashed lines 
are linear fits for all bond dimensions D for sizes L — 8, 16 
respectively. 
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FIG. 7: Overlap of TPS of various bond dimensions D 
with the exact ground state wavefunction |^>)o obtained by 
exact diagonalization for system size L = 4. For D — 16 the 
overlap is 0.99979. 
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L = 16, the z component is zero for all available bond 
dimension D. Asymptotically, as D increase, one could 
expect C z (L/2,L/2) grows to different values for differ- 
ent system sizes L, and for sufficiently large system size 
C z (L/2, L/2) —> due to automatic symmetry breaking. 

In Fig. [7] we calculate the overlap of the TPS wavefunc- 
tion for various bond dimensions D with the exact 
ground state | -0) o obtained by an exact diagonalization 
for a 4 x 4 system. For D = 16 the overlap is 0.99979. 

For all the data presented, the maximum number of 
singular values kept at each iteration step during the con- 
traction is D cut = 2D for all bond dimensions D. 



FIG. 5: Absolute error of the energy per bond as a function 
of the bond dimension D for system sizes L = 4, 8, 16 on a 
normal scale. 

























A" 


+ A L * 




CvJ 


0.05 






A - -Jt" 


L=4 C x ' y 




_i 










L=8 C x ' y 




cm" 
o 




*' 






L=16C x ' y 
L=4 C z 
L=8 C z 


i— Ah 






* 




. . ' 


L=16C Z 
. • • * 


» - • 




0.00 


















5 


10 




15 



D 

FIG. 6: Spin-spin correlation C a (L/2, L/2) (a = x,y,z) as 
a function of D for system sizes L — 4, 8, 16. Solid lines are 
the x, y components, dash lines show the z component. Note 
that C x = C v for the [7(1) symmetry. 



IV. SUMMARY AND DISCUSSION 

In this paper, we proposed a variational Quantum 
Monte Carlo (vQMC) algorithm to evaluate a tensor 
network state of relative large bond dimensions. We il- 
lustrated the Monte Carlo (MC) sampling procedure in 
terms of the tensor entanglement renormalization group 
(TERG) contraction algorithm. We applied this method 
to the well studied Antiferromagnetic Heisenberg model 
on square lattice. Upon obtaining the ground state wave- 
function via imaginary time evolution with essentially no 
notion of the lattice size, we evaluated the ground state 
energy and the staggered magnetization for systems on 
finite square lattices using MC sampling method. Not 
surprisingly, we found that the converged tensors ob- 
tained for the infinite lattice give very accurate results 
also when used for considered finite size lattices. The 
wavefunction of a finite bond dimension D thus can be 
used to reliably extrapolate the expectation values in the 
thermodynamic limit through finite D scaling followed 
by finite size scaling. We have shown that the tensor 
network ansatz based vQMC method is a promising way 
to go to a very large bond dimension and thus allowing 
reliable study of many interesting models. 
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We do not claim that the tensor product state (TPS) 
describing the ground state wavefunction for an infinite 
lattice obtained via the simple update is the ultimate 
solution for a finite lattice. One still need further op- 
timization for a finite lattice if initializing from a TPS 
describing the infinite lattice. Many previous studies 
of 1 dimensional systems had used the matrix product 
state (MPS) obtained from an infinite chain algorithm 37 
to initialize the optimization for a finite chain and ob- 
tained remarkably good result s 38 ' 39 . We will discuss how 
to update tensors for a finite 2 dimensional system else- 
where^. Another advantage of the MC sampling method 
is the possibility of incorporate lattice and spin symme- 
tries into the MC sampling scheme to improve accuracy, 
which has been demonstrated in Ref4^ in the case of 
simulating 2 dimensional system via scale-renormalized 



MPSs. Since we obtained tensor directly from the simple 
update, we did not employ symmetries in the sampling 
procedure. 
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